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Abstract 



J^ , The two-dimensional flow of viscous incompressible liquid in a square cavity with 

Q ■ a free boundary and differentially heated vertical sides is considered in the present 

Q ! work. The influence of gravitational and thermocapillary convection on temperature 



and velocity fields is studied in large range of dimensionless parameters and similarity 



^ , criteria using equations in a Boussinesq approximation. Limiting cases of dimension- 

c/3 ' less parameters are analyzed numerically. 

Ph. 

_ Introduction 

X 

c^ I It is known that in the nonuniformly heated fluid motion appears. Without free 

boundaries it appears due to thermal (gravitational) convection. In the nonuniformly 
heated liquid with a free boundary there arises thermocapillary convection |[T|, |0. 
Under the low gravity condition it is necessary to take into account both of these 
processes. It is interesting to analyze the influence of these factors on heat fluxes 
and velocity fields formation under the non-gravity condition and with an increase in 
gravity. 

We will consider a plane stationary convectional system with a side heating. Such 
a formulation of the problem is realized in practice, for example, in crystal growth 



equipment, various energy plants [0]) 0- That is why the knowledge of flow and heat 
fluxes structure in the large dimensionless parameters range (Rayleigh Ra, Prandtl 
Pr and Marangoni Ma numbers) is of scientific and practical interest. 

Model 

The stream function - vorticity {ip—uj) formulation of the problem was used. For plane 
geometry initial non-dimensional stationary convection equations in a Boussinesq 
approximation under the non-gravity condition are the following [Q: 

-M + -M = ^v-, (1) 

V¥ = ^, (2) 

^{uT) + ^{vT) = -^V^T. (3) 

dx^ ^ dy^ ^ Ma ^ ^ ^ 

where « = |^, v = -^, u = ^ - ^. 

ay ' ax' ay ox 

The problem is characterized by the following parameters: the Marangoni number 
Ma = — ^^— and the Prandtl number Pr = —. Here AT = Th — Trj, Tu is the 

/xa pa n K^ 1 n 

hot wall temperature, Tc is the cold wall temperature, ar denotes the temperature 
coefficient of surface tension, L stands for the characteristic length (the side of a 
square cavity), /i is the dynamic viscosity, p designates the density, a is the thermal 
diffusivity. 

Calculations with various Prandtl and Marangoni numbers were carried out. Two 
limiting cases were also considered: Ma -^ and Pr -^ oo. li the Marangoni number 
tends to zero (i.e. surface tension force is equal to zero) then the set of Eqs. (|l| - ^ 
is modified as follows: 

V'^ = 0, (4) 

V'V' = ^, (5) 

V^T = 0. (6) 

If the Prandtl number tends to infinity (that is the case of the strongly viscous fluid) 
then the set of Eqs. (|I| - 0) is developed as follows: 

V'^ = 0, (7) 



V''* = w, (8) 

If the force of gravity g is not equal to 0, initial stationary non-dimensional equa- 
tions of thermal convection in a Boussinesq approximation in the uniform gravita- 
tional field in [ip — uj) variables are then |Q: 

Pr o 
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(12) 



Here Ra = — — ^^"~ '^' is the Rayleigh number, (3 = —-^ denotes the volumetric 
coefficient of thermal expansion. 

In the present work the following boundary conditions are considered. The vertical 
sides are at temperatures T = Th = 0,5 ("hot" wall) and T = Tc = —0,5 ("cold" 
wall). The lower horizontal wall and the free surface y = 1 are insulated. Both 
velocity components are zero on the walls. On the free surface the component /i|^ 
of the viscous tension tensor is equal to the tangential force acting on the surface 
—ct^, and vertical velocity v is equal to zero. Thus the boundary conditions for the 
considered equations can be written in the following way: 

0, T = Th = 0,5 at X = 0, 

0, T = Tc = —0,5 at X = 1, 

0, f = at y = 0, 

^ = 0^ = -f> f = o -t, = i. 

The case of the surface tension force equal to zero was also considered. This 
corresponds to the system (|5 - ||) with the boundary condition |^|j/=i = : 



^ = 


dtp 
" dx 


iJ = 


" dx 


iJ = 


dtp 
■ dy 



^ = 0, = at y = 1. 



Method of solution 



Formulation of the problem (|l] - §), as many other problems of viscous incompressible 
fluid in {ip,uj) variables has the following difficulty. Boundary conditions on the walls 
are determined only for the stream function, and not for the vorticity, which is defined 



only inside the domain according to the Eq. (^. To overcome this difficuhy various 
approaches are used, for example, approximate boundary conditions for vorticity. In 
the present work we use the Toma condition for uj on the wall 0: 

^^^2(*^L_1J«+0(A;,). (13) 

Here A/i denotes the mesh size, tpk is the value of the stream function in the boundary 
node fc, ipk+i is is the value of ip in the node k + 1 nearest to the wall. 

Calculations were also performed with the second order boundary condition for 
uo, namely. Woods condition. Computational results with these different conditions 
are almost the same. However, the use of the approximate boundary conditions for 
vorticity on the wall for the Eqs. ([^ - ^ at high Marangoni numbers and on fine 
meshes leads to the considerable slowing-down of the convergence. That is why we 
also used the boundary conditions calculation method @], which allowed us to improve 
the convergence in 2 — 4 times, and in some cases on the order in comparison with 
the use of the Toma formula. 

Idea of this method is to determine the boundary condition for vorticity inside 
the main domain, where u is defined according to (0). Equation for uj (||) is solved in 
the auxiliary domain. The solid boundaries of this domain are displaced on the mesh 
size into the cavity from the solid boundaries of the main domain. The free surface 
y = 1 is common for these domains. On the free surface vorticity is determined as 
follows: 

, ,1 — du\ _ _&T\ 

Stream function and temperature equations are solved in the main domain. There 

are two boundary conditions for the stream function. The condition ^/'o = is used to 

solve equation for ip. The obtained stream function field does not satisfy the gradient 

condition (^)|2,=o,a=i,j/=o = yet. Therefore, the values of ip on the boundary of 

the auxiliary domain are corrected with the help of the difference analog of this 

condition |^. Using a three-point approximation of the second order for the derivative 

(|^)U=o,a;=i,j/=o we will obtaiu: 

at a; = ^2j = ^i'sj, J = 2, . . . , m - 1; 

at a; = 1 i^n-ij = \ipn~2j, j = 2, . . . , m - 1; 

at y = ^i2 = 2^i3, i = 2, . . . , n - 1. 



When approximating Eqs. (|1|, H) and (|T^, |T^ for temperature and vorticity ex- 
ponential fitting discretization (or scheme of integral identities) was used 0, 0. It 
allowed us to obtain a higher precision in comparison with the usual approxima- 
tions. As a result, a five-point algebraic system was carried out. Equations for uj 
and T do not satisfy the diagonal dominance condition. It is known that without 
this condition, many effective methods fail to converge or converge very slowly. In 
the present work modification of Buleev's method §] and a splitting method [jlO| are 



used. They enable to find a solution of the system without a diagonal dominance 
property. Buleev's method modification converges faster than the splitting method 
when the stream function and the temperature equations are solved. But for the 
vorticity equation with high Marangoni numbers {Ma > 10^) Buleev's method does 
not converge to the necessary precision (namely, l.e — 7). Thus the equation for u 
was solved by splitting which enables to obtain the prescribed precision. 

To improve convergence the damping operation for vorticity was used. It is de- 
termined by the following recurrent relation: 

uT' = Ouj-^' + (1 - 0)u;,^ 

where 9 stands for the damping parameter, w^" is the damped value of vorticity from 
the n iteration, o;""*"^ denotes the value of u from the n + 1 iteration. 

During numerical experiments for a 61 x 61 mesh the damping parameter 9 was 
about 0,002. 

Main results and their analysis 

We use rectangular non-uniform thickening to the boundary of the 21 x 21, 41 x 41 
and 61 x 61 area grids. Calculations were performed at the Marangoni numbers from 
10^'^ to 10^, the Prandtl numbers from 1 to 100 and the Rayleigh numbers from 
to 10^. Two limiting cases: Ma -^ and Pr -^ oo were also considered. Results 
are presented in Figures (1-10) and in Tables I and II on the 61 x 61 mesh. Results 
at the Marangoni number Ma = 10^ are not very reliable as the calculating scheme 
converges badly at this value of dimensionless parameter. 

Figure 1 shows the infiuence of surface tension force on the temperature field 
under the non-gravity condition {Ra = 0). If the Marangoni number Ma = 0, i.e. 



the surface tension force is equal to zero (Fig. la), there is no convection, and the fluid 
rests. With an increase in the Marangoni number the convective mixing enhances, 
and at Ma = 10^ the flow becomes vortical. One can see from Figures lb, Ic, Id, 
that the greater Ma, the stronger the contour maps deflect and are pressed to the hot 
and cold walls. Near these walls temperature boundary layers appear. The boundary 
layer near the hot wall is wider than near the cold. One characteristic property of the 
picture is the condensation of contour maps to the upper right corner of the cavity, 
to the region near the free surface and the cold wall. In the center of the domain 
a vortex is formed. Contour maps on Figure 1 correspond to the Prandtl number 
Pr = 1. Figure 2 shows temperature contour maps at Pr = oo (that is the case of the 
strongly viscous fluid, Eqs. (0 - §))• Behavior near the free surface y = 1 is different 
for the corresponding contour maps of Figures 1 and 2. In the case of the strongly 
viscous fluid contour maps come normally to the free surface and remain undeflected 
near it longer than in the case of Pr = 1. The explanation is that it is more difficult for 
the surface tension force to mix more viscous fluid. From the comparison of contour 
maps at Ma = 10^ (Figures Id and 2c) one can see that at Pr = 1 a vortex region is 
larger and extend almost over the entire domain. At Pr = oo the vortex is situated 
closer to the free surface, convective flow in the entire domain is weaker. Thus the 
more viscosity of the fluid, the weaker convection there at the same surface tension 
force. 

Figure 3 presents the influence of the gravity (Ra) on the temperature field. An- 
alyzing the results, we can say that with the increase in Rayleigh number the whole 
picture of the flow becomes more complex, convection mixing of the fluid intensifies, 
the vortical structure of the stream and narrow boundary layers with the sharp tem- 
perature difference across the boundary layer near the hot and cold walls appear. 
Comparison of corresponding contour maps of Figures 1 and 3 shows that the surface 
tension force tends to form a vortex in the center of the domain and to condense the 
contour maps to the upper part of the cold wall; whereas the force of gravity tends 
to form a vortex closer to the boundary of the cavity, remaining its center without 
a vortex, and to condense the contour maps to the vertical walls. At Ra = 10^ two 
narrow boundary layers near the vertical sides of the width ^ 0,07 each, with hori- 
zontal temperature gradient are generated. In the remaining part of the calculation 
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102 


11,871 


0,6776 
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103 


13,639 


0,6776 
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10^ 


26,626 


0,5 


1 


10^ 


73,827 


0,3225 


1 


10^ 


192,88 


0,0541 


0,5969 



Table 1: Dependence of the maximum value of velocity module and its location on 
the Rayleigh number. 

domain of the width ^ 0,85 the contour maps are located at a greater distance and the 
temperature gradient points vertically upwards. Such behavior of the temperature 
contour maps influences the stream function field. One can see from Figure 4d for the 
ip contour maps {Ra = 10^) that near the vertical walls contour maps condense, and 
two narrow boundary layers with the sharp difference of stream function values across 
the boundary layer are formed. In the remaining part of the calculation domain the 
contour maps are located at a greater distance. At Ra = 10^ (Fig. 4c) the secondary 
flows are formed in the center of the domain. 



Table I shows an increase in the maximum value of velocity module \^ii? + v^ and 
the change in its location with an increase in the Rayleigh number. At Ra < 10^ this 
maximum is situated on the free surface, and the greater Ra, the closer the maximum 



to the hot wall. At Ra = 10^ the maximum value of ^/v? + v^ is located close to the 
center of the hot wall. 

Dependence of the maximum value location of ^/u"^ + v"^ on the Marangoni number 
was built (Fig. 5). For all considered Ma numbers the maximum value of velocity 
module is situated on the free surface. With an increase in Ma it moves to the upper 
right corner, the common point of the cold wall and the free surface. For a more 
viscous fluid this maximum is located closer to the hot wall for the same Marangoni 
numbers. It is evident that the surface tension force tends to move the vortex to 
the upper boundary of the cold wall, and it is more difficult to mix more viscous 
fluid. In the range of Ma from 10^ to 10^ some oscillation of the curve can be seen 
(Fig. 5). It is connected with the more complex structure of the flow at the high 



Marangoni numbers. A high precision is necessary for the calculation of the flow at 
the Marangoni numbers Ma > 10^. Such precision could not be obtained with the 
help of the proposed scheme. 

Figure 6 presents the contour maps of stream function at Pr = 1, Ra = for 
various Marangoni numbers. With an increase in Ma the contour maps begin to 
come closer to the upper corners, and in the lower corners the secondary flows appear. 
With the increase in the Marangoni number the inner vortex domain increases, and 
the contour maps value module corresponding to this vortex decreases. The contour 
maps condense near the free surface. 

Figure 7 shows ip contour maps for the strongly viscous fluid at different Marangoni 
numbers. One can see from a comparison of Figures 6 and 7 that in more viscous 
fluid the convective flow in the entire calculation domain is weaker, vortex is located 
closer to the free surface and vortex's region is less; there are no secondary flows in 
the lower corners of the domain. 

Let us compare Figures 6a {Pr = 1, Ma = 10^, Ra = 0) and 4a {Pr = 1, Ma = 
10^, Ra = 10^). With the appearance of the thermal gravitational convection the 
contour maps shape changes. The lower parts of the contour maps are located closer 
to the bottom of the domain, and contour maps tend to the lower corners. It is the 
result of the gravity force influence. The greater Ra, the stronger the influence of 
the gravity force in comparison with the surface tension force on the character of the 
flow. Thus at Ra = 10^, Ma = 10^ (Fig. 4) the contour maps also tend to the upper 
right corner, i.e. the thermocapillary convection considerably influences the contour 
maps shape. With the further increase in Ra the gravitational convection becomes 
dominant. Beginning from Ra = 10^ the secondary flows are formed. One of them 
tends to the upper left corner of the cavity, to the hot wall, another - to the lower 
right corner, to the cold wall (Fig.4b, 4c, 4d). 

Calculations of the convection problem under the non-gravity condition with the 
surface tension force equal to zero (Eqs. (^ - ^) with the zero boundary condition 
^\y=i = 0) were made. Here no forces operate on the fluid, there is no convection, 
and all profiles of stream function, velocities and vorticity at various sections of the 
calculation domain are almost zero. The temperature contour maps in this case are 
straight lines (Fig. la). It means that the heat is transferred only by heat conduction. 



Calculations of the convection problem without gravitation were also carried out 
with the following parameters: Ma = 1; 0,1; 0,01; 0,001; Pr = 1. With the decrease 
in the Marangoni number the stream function and velocity profiles tend monotonously 
to the zero values, to the profiles of the limiting case Ma = 0. The temperature 
contour maps also tend monotonously to the contour maps of the case Ma = with 
the decrease in Ma. 

Another limiting case: Pr -^ oo for various Marangoni numbers (Eqs. (0 - |^)) 
under the non-gravity condition was also considered. Computational results show 
that for every Marangoni number between and 10 the flow does not depend on the 
Prandtl number (at Pr between 1 and cxo). At Ma = 10^ and Pr from 10 to cxd the 
corresponding profiles of the flow agree (Fig.2a (temperature contour maps). Fig. 7a 
(contour maps of ^)), but differ from the flow profiles at the Prandtl number Pr = 1 
(Fig. lb, 6a). A similar situation takes place at Ma = 10'^ (Fig.lc, 6b; 2b, 7b). At 
Ma = 10^ and Pr = 100, oo the profiles are close to each other and differ from the 
flow profiles at Pr = 10. 

For the heat balance control the average Nusselt numbers on the hot and cold 
walls and on the vertical mid-plane of the cavity were calculated. The local Nusselt 
number Nu{x,y) at the point {x,y) and the average Nusselt number Nux^ on the 
section x = Xq are determined as follows: 

dT 

Nu{x,y) = Ma-uT — -—, 

LJJb 

Nu^o = / Nu{xo,y)dy. 
Jo 

Table II shows the computational results of the convection problem in the cavity 
with a free surface and side heating under the non-gravity condition at Pr = 1 and 
various Marangoni numbers. Here ^(1/2, 1) is the value of m velocity at the mid-point 
of the free surface, Nuq, Nui/2, Nui designate the average Nuselt numbers on the 
sections x = 0, x = l/2, x = 1 respectively. The obtained results were compared 
with [|lT]. In the work |jTl| a high accuracy scheme was used to solve Eqs. (|l] - ^. 



One can see from the Table that the present work method allows us to obtain with 
acceptable precision the integral characteristics of heat-removing {Nuo, Nui/2, Nui). 
However, the local characteristics {Nu{0, 1), Nu{l, 1)) are determined with essential 



Ma 


grid 


u(l/2,l) 


Nuq 


Nui/2 


Nui 


Nu{^,l) 


Nu{l,l) 


W 


21 X 21 
41 X 41 
61 X 61 


1,0326(-1) 
1,0711(-1) 
1,0671(-1) 


1,1322 
1,1010 
1,0945 


1,0999 
1,0977 
1,0953 


1,0396 
1,0867 
1,0924 


0,7972 
0,7402 
0,7313 


1,7816 
1,9874 
1,9905 


[11] 


61 X 61 


1,0869(-1) 


1,0962 


1,0962 


1,0962 


0,7301 


2,028 


103 


21 X 21 
41 X 41 
61 X 61 


5,2221(-2) 
5,2032(-2) 
5,1404(-2) 


2,1800 
1,9858 
1,9437 


1,9976 
1,9742 
1,9528 


1,6900 
1,8983 
1,9297 


1,3444 
1,0268 
0,9707 


7,3224 
10,6884 
11,438 


[11] 


61 X 61 


5,0018(-2) 


1,9258 


1,9258 


1,9258 


0,9550 


11,75 


10^ 


21 X 21 
41 x41 
61 X 61 


2,7263(-2) 
3,1541(-2) 
3,1562(-2) 


5,4250 

4,5622 
4,4459 


3,9362 

4,3744 
4,4467 


2,6086 
4,1286 
4,3524 


5,1195 
3,1227 
2,5091 


12,6693 
38,5290 
60,471 


[11] 


61 X 61 


3,0381(-2) 


4,3621 


4,3621 


4,3654 


2,2334 


77,09 



Table 2: Computational results at Pr=l. 



errors, particularly at Ma = 10^. 

Dependencies of the average Nusselt numbers on the Prandtl, Marangoni and 
Rayleigh numbers (Fig. 8, 9 and 10 respectively) were also calculated. Fig. 8 shows 
that the greater the Marangoni number for the same Prandtl number (i.e. the 
stronger thermocapillary convection), the stronger heat-removing. At Ma between 
and 10^ the heat-removing coefficient does not change with the increase in Pr. At 
Ma = 10^; 10^ and from Pr = 10 this coefficient takes some constant value, and 
every Marangoni number has its own value. Besides, at Ma between 10^ and 10*^ an 
abrupt increase in the Nusselt number is observed, especially for the case Pr = 1. 
Probably, the reason for it lies in the appearance of a more complex structure of the 
flow and in the loss of it the stability [^ , which lead to the sharp increase in heat- 
removing. Such an abrupt increase in heat-removing is also observed in Fig.9, which 
presents dependence of Nu on Ma for various Prandtl numbers. The Nusselt number 
increases with the growth of the Rayleigh number, which characterizes the intensity 
of the gravitational convection (Fig. 10). At Ra between 10^ and 10^ abrupt increase 



in the Nusselt number also occurs because of the appearance of a more complex flow 
structure. Figures 9, 10 present functions that are the parameter approximations of 
dependences of the Nusselt number on the hot wall on the Rayleigh number at Ra 
between 10^ and 10^ and on the Marangoni number at Ma between 10^ and 10^. In 
figures these functions are showed by triangles (AAA). They are the linear combi- 
nations of functions with the fitting coefficients as follows: 

g{Ma) = 1,481 + 7,15 ■ 10"^ • Ma - 3,277 ■ 10"^ ■ Ma"^ - in the case of the Nusselt 
number on the Marangoni number dependence; 

g{Ra) = 4,122 + 9,507 • 10"^ ■ Ra - 4,296 ■ lO^^^ ■ Ra"^ - in the case of the Nusselt 
number on the Rayleigh number dependence. 



Conclusions 

1. The influence of gravitational and thermocapillary convection on the temperature 
and velocity fields has been studied. Computational results show that at Ma < 10^ 
and beginning from Ra = 10^ the gravitational convection becomes dominant. 

2. With the increase in the gravity force the whole picture of the flow becomes more 
complex, convection mixing of the fluid intensifles, narrow boundary layers with the 
sharp temperature and velocity difference across the boundary layer close to the hot 
and cold walls appear. Beginning from Ra = 10^ the secondary flows formation 
occurs. 

3. In a more viscous fluid the convective flow in the entire calculation domain is 
weaker, vortexes are located closer to the free surface and vortexes region is less. 

4. Two limiting cases: Ma — > and Pr -^ oo without gravitation were considered. 
Computational results show that with the decrease in the Marangoni number the flow 
profiles tend monotonously to the flow proflles of the case Ma = 0. Ranges of the 
dimensionless parameters where the flow does not change were also determined. 

5. The parameter approximations of dependences of the Nusselt number on the hot 
wall on the Rayleigh number at Ra between 10^ and 10^ and on the Marangoni 
number at Ma between 10^ and 10^ were built. 
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FIGURES 

Fig.l. Contour maps of temperature T : Pr = 1, Ra = 0, 

a - Ma = 0, b - Ma = 10^ c - Ma = 10^ d - Ma = 10^. 

Fig. 2. Contour maps of temperature T : Pr = oo, Ra = 0, 

a - Ma = 10^, b - Ma = 10^, c - Ma = 10^ 

Fig. 3. Contour maps of temperature T : Pr = 1, Ma = 10^, 

a- Ra = 10^, b - i?a = 10^ c - Ra = 10^ d - i?a = 10^. 

Fig. 4. Contour maps of stream function ip : Pr = 1, Ma = 10^, 

a- Ra = 10^, h- Ra = 10^ c - Ra = 10^ d - i?a = 10^. 

Fig. 5. Maximum value location of the velocity module versus the Marangoni number 

(— Pr = 1, -- Pr = cx)). 

Fig. 6. Contour maps of stream function ip : Pr = 1, Ra = 0, 

a - Ma = 10^, b - Ma = 10^ c - Ma = 10^ 

Fig. 7. Contour maps of stream function ip : Pr = oo, Ra = 0, 

a - Ma = 10^, b - Ma = 10^ c - Ma = 10^ 

Fig. 8. Average Nusselt number on the hot wall versus the Prandtl number 

(1 - Ma = 0, 2 - Ma = 10^, 3 - Ma = 10^ 4 - Ma = 10^). 

Fig. 9. Average Nusselt number on the hot wall versus the Marangoni number 

{g{Ma) = 1,481 + 7,15 ■ 10"^ ■ Ma - 3,277 ■ 10"^ ■ Ma^) 

— Pr = 1, Pr = CX), A A g{Ma). 

Fig. 10. Average Nusselt numbers on the hot and cold walls versus the Rayleigh 
number {g{Ra) = 4,122 + 9,507 ■ 10"^ ■ Ra - 4,296 ■ lO^^^ . ^^2^ 

— iVwO, --iVwl, AA g[Ra). 



o 



CM 
OO 



CO 
o 



CO 
CO 



CO 

eg 
eg 
o 



in 
in 

o 



CD 

CO 



o 

-;^ 

CJl 
CJl 



CO 
t^ 
rg 
rg 



(J) 
o 



rg 

00 



CO 
o 



/ /// 


/// 


7 


f / /// 


/ / 


/ / / 


/ -" / 


' 


1 / 9 / / / 


^ s? / 


1 / 




/ / # / 


r / 


eu / 
- / 




^ / 


/ / 


« 




/ o / 


/ / 


9 




1 

eg 


/ / 


/ / 


/ 




CO 




1 lO 


/ / 


/ 




T 




in 


/ / 


/ 




00 




V 


/ 


1 




o 




o 


Cv^ 












o 


r^ 










o 




eg 
















eg 










CO 






<D 


o 


















'^ 












































I 


Y^ -^y^/ll 


C 

c 






o 


I .0--'^/// 








f /'"'^// 




p 




/ <^ 




" , / F / 




i ^ 




\ \ 


■* 




\\\ \ 


o 



b 



/J 




/ yjM 


^0-7364^ 


^ y — 


///n\ 




V 


^ 


f 


o 
to 

tvj 

CO 

1 


.0.1364 . 




11 




^^■^" / 


f o 

/' 


r 


1 






10° 10' 10' io\jo^ 

Ma 



